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Abstract. Hydrodynamic (Landau) instability in combustion is typically ass^.a.^^thtte 

rate and tlJ disturbance wavenumber, significantly generalizing previous classical results for this 
m-oblem thaUissumed a constant norma, burning rate. Additionally, however a ° 

thresholds * It is found that for sufficiently small values of this parameter, there exists a stable 
the pulsating stability boundary evolves into a C-shaped curve m 

oressure-sensitivity) plane, indicating loss of stability to pulsating perturbations for all sutticenl y 
We IS wavelengths. It is thus concluded, based on charactenstic parameter values that 
Tn eaX I form of hydrodynamic instability in liquid-propellan, combustion is of a non- 

eZ n w vVn tore, distinct from the steady, cellular form originally predicted by Landau. 


1. Introduction 

The notion of hydrodynamic instability in combustion originated with Landau’s seminal study 
propagation (Landau ID). In .ha, work. „ was postulated rha, a flam. could 
be represented by Lmface of discontinuity, with an appropriate density jump across the flame 
surface, and the normal equations of hydrodynamics governing the flow on either sirk ^of th 
flame surface Assuming that this surface propagated normal to itself with constant speed, 
was then determined through a straightforward linear stability analysis that a P« ^gaseous 
flame was intrinsically unstable to steady (cellular) disturbances wheneve as s g^emlly *e 
case the density on the burnt side of the front was less than that of the unbumt mixture. Smce 
this result is at variance with experimental observations of steady, P ' a ^ 
has sparked numerous other theoretical investigations that have sought to resolve this paradox, 

* This work was supported by the US Department of Energy under contract DE-AC04-94AL85000 and by the NASA 
Microgravity Science Research Programme under contract C-3203 l-fc. 
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through both the inclusion of additional effects and the analysis of less phenomenological 
combustion models. Today, this specific form of hydrodynamic instability, generally referred 
to as Landau instability, plays a central role in the study of flame/flow interactions. Though 
Landau’s analysis essentially collapsed all thermal effects to the change in density across the 
front, Landau instability is nonetheless implicitly associated with the thermal expansion of the 
gas arising from the heat release during combustion. Indeed, this instability persists even when 
the entire flame structure is resolved (cf Matalon and Matkowsky [2], Pelce and Clavin [3]), 
although it is now recognized as a predominantly long-wave cellular instability that is subject 
to various stabilizing influences for disturbances corresponding to larger wavenumbers. 

A related problem briefly addressed at the end of Landau’s original study (for which a 
more standard and less controversial type of result was obtained), and the one of interest here, 
was a modest modification of the flame-propagation problem just described. In this modified 
problem, the unbumt mixture was taken to be a liquid propellant while the burnt region again 
consisted of gaseous products. The physical existence of a liquid/gas interface led to the 
inclusion of additional physics in the model, namely surface tension at the interface and the 
influence of gravitational acceleration (downward propagation was assumed). As a result, and 
in contrast to strictly gaseous flame propagation, a stability criterion was derived such that the 
liquid/gas interface was either hydrodynamically stable or unstable in the Landau (cellular) 
sense depending on whether the product of the gravitational acceleration and the coefficient of 
surface tension was greater or less than a critical value. This result was later extended by Levich 
[4], who considered the effects of (liquid) viscosity in lieu of surface tension and obtained a 
similar result for the product of the gravitational acceleration and the viscosity of the liquid. 
These classical models of liquid-propellant combustion, despite their relative simplicity, have 
long served as seminal examples that correctly describe the onset of (cellular) hydrodynamic 
instability in reactive systems. 

Although the assumption of a thin reaction region remains a frequently valid and 
useful simplification in many combustion studies, the classical assumption of a constant 
normal burning rate has long been regarded as a probable oversimplification when applied 
to the problem of combustion instability. Early attempts at modification began with the 
phenomenological assumption of a linear relationship between the burning rate and the flame 
curvature (cf Markstein [5]), while more modem approaches have employed sophisticated 
asymptotic methods to analyse the flame structure and to derive formal expressions for the 
burning rate as a function of local conditions at the combustion front, for both gaseous 
and condensed combustion waves (cf [2,3], Sivashinsky [6], Margolis [7]). In propellant 
combustion, on the other hand, it has long been customary to experimentally measure the 
pressure response, or pressure sensitivity, of the burning rate, as well as (to a lesser extent), 
its temperature sensitivity. Although asymptotic models that resolve the combustion-wave 
structure can be developed, under various approximations, for propellant combustion (cf 
Margolis and Williams [8, 9], Bechtold and Margolis [10]), it nonetheless remains true that the 
direct assumption of a combustion surface that propagates according to a prescribed burning- 
rate law, such as a phenomenological law based on experimental observations, allows one to 
circumvent the intricacies of the combustion region and to impose fewer restrictive assumptions 
on the hydrodynamic model for the unbumt and burnt regions. 

Accordingly, for the case of liquid-propellant combustion, the classical Landau/Levich 
hydrodynamic models may be combined and extended to account for a dynamic dependence, 
absent in the original formulations, of the instantaneous burning rate on the local pressure 
and temperature fields (Armstrong and Margolis [II, 12]). Although exact expressions for the 
neutral stability boundaries in question can be obtained in specific limiting cases, more general 
results may be obtained by exploiting the realistic smallness of the gas-to-liquid density ratio 
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p. Thus, neglecting thermal coupling effects, an asymptotic analytical expression was derived 
for the cellular stability boundary A„(k), where A p is the pressure sensitivity of the burning 
rate and k is the wavenumber of the disturbance (Margohs [13]). The results demonstrate 
explicitly the stabilizing effect of gravity on long wave disturbances, the stabilizing effec 
of viscosity (both liquid and gas) and surface tension on short-wave perturbations and he 
instability associated with intermediate wavenumbers for critical negative values of A p . In the 
limiting case of weak gravity, it was shown that hydrodynamic instability in liquid-propellant 
combustion is a long-wave instability phenomenon, whereas at normal gravity, this Jnstabdity 
is first manifested through 0(1) wavenumber disturbances. It was also demonstrated that, 
in general, surface tension and the viscosity of both the liquid and gas phases each produce 
comparable stabilizing effects in the large- wavenumber regime, thereby providing impo an 
modifications to the previous classical analyses in which one or more of these effects was 
neglected For A p = 0, the Landau/Levich results are recovered from the new model in 
appropriate limiting cases, although, in the realistic limit of small p. this corresponds to a 

hydrodynamically unstable parameter regime [1,4, 13]. 

In addition to the classical cellular form of hydrodynamic stability described thus far, there 
exists in this extended model a pulsating form of hydrodynamic instability corresponding to the 
loss of stability of the steady, planar solution to time-dependent perturbations (Margohs [14]). 
This occurs for negative values of the pressure-sensitivity parameter A p , and is thus absent 
from the original Landau/Levich models. Nonetheless, its existence in the extended model 
results in a stable band of negative pressure sensitivities bounded above by the classical Landau 
type of hydrodynamic instability, and below by the pulsating form of hydrodynamic instability 
just described. A stable range of negative pressure sensitivities is applicable, for example, to 
certain types of hydroxylammonium nitrate (HAN)-based liquid propellants at low pressures 
for which non-steady modes of combustion have been observed, but not otherwise investigated 
further (Vosen [15]). While such non-steady combustion may correspond to secondary and 
higher-order bifurcations above the cellular boundary (cf Bechtold and Margohs [16]), it may 
also be a manifestation of the pulsating form of hydrodynamic instability just described. 

In the present work, we explore the effects of incorporating a non-zero temperature 
sensitivity into our previous asymptotic analyses in the limit of small gas-to-liquid density 
ratios [13, 14}. This entails a coupling of the energy equation for temperature to the previous 
purely hydrodynamic problem, and leads to a significant modification to the pulsating boundary. 
Specifically, for sufficiently large values of the temperature- sensitivity parameter, this boundary 
evolves into a C-shaped curve such that liquid-propellant combustion can become intrinsically 
unstable to this newer form of hydrodynamic instability. 

2. The mathematical model and its basic solution 

The mathematical model was described previously [13], but is briefly summarized here for 
completeness. In particular, it is assumed, as in the classical models, that there is no distributed 
reaction in either the liquid or gas phases, but that there exists either a pyrolysis reaction or an 
exothermic decomposition at the liquid/gas interface that depends on local conditions there. 
For simplicity, it is assumed that within the liquid and gas phases separately, the various 
physical properties are constants, with appropriate jumps in these quantities across the phase 

boundary. . . 

The non-dimensional location of the downward-propagating hquid/gas interface as a 
function of space and time is denoted by x-s — d 2 * * 5 v (x | , x 2 < 0. where xi is the vertical coordinate 
and the adopted coordinate system is fixed with respect to the stationary liquid at xj = 

Then, in the moving coordinate system x = x u y = x 2 , z = x 3 - <t>, (xi,x 2 ,r), in terms o 
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which the liquid/gas interface always lies at z — 0, the complete formulation of the problem 
is given by conservation of mass, energy and momentum within each phase as 


V ■ u — 0 z/ 0 


( 1 ) 
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where v , 0 and p denote velocity, temperature and pressure, respectively, Pn and Pr g denote 
the liquid and gas-phase Prandtl numbers, p, A. and c (used below) are the gas-to-liquid 
density, thermal diffusivity and heat-capacity ratios, and Fr is the Froude number. These non- 
dimensional variables have been defined in terms of their dimensional counterparts (denoted 
by tildes) as 
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where U is the reference propagation speed of the interface for the case of steady, planar 
deflagration. The non-dimensional parameters, some of which first appear below in the 
conditions at the gas/liquid interface, have been defined as 


P = 


h 

pi 





(5) 


where T u is the unbumt (liquid) temperature at z = — oo and T a is the adiabatic burnt (gas) 
temperature at z = +oo, y is the coefficient of surface tension for the liquid surface, and \>i and 
v g are the liquid and gas-phase kinematic viscosities. We note that the inverse Froude number 
Fr~ l represents the non-dimensional gravitational acceleration and that pXPr g = piPri , 
where p = p g /p.i is the gas-to-liquid viscosity ratio. The non-dimensional mass burning rate 
appearing in equation (9) below is defined by A(p\ z =o> y 0|*=o) = T\ h= ^)/p { U y 

and is assumed to be functionally dependent on the local pressure and temperature at the 
liquid/gas interface. By definition, A = 1 for the case of steady, planar burning, but 
perturbations in pressure and/or temperature result in corresponding perturbations in the local 
mass burning rate. 

Equations (l)-( 3) are subject to the boundary conditions 


v = 0 



at z = —oo 
at z — +oo 


0|^=o- = 0lz-o + 


(6) 


and appropriate jump and continuity conditions at the liquid/gas interface. The latter consist 
of continuity of the transverse velocity components (no-slip), 


n s x = n s x v 


(7) 
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where v± = v| I=0 ±, conservation of (normal) mass flux, 

3<t>, 

h s • (u_ — pv+) = (1 — p)5(4>j) — - 

the mass burning rate (pyrolysis) law, 

n 5 • v- - SOD,)^ 1 = Mp\z=o>, ©U=o) 
ot 

conservation of flux of the normal and transverse components of momentum, 
plz=o- - p\tm 0* = n s • [pv + (n s ■ V + ) - »-(«, • v .) - pkPr g e+ ■ n s + Pr,e _ • n s ] 

r 9<i>, 1 

n, x pu + (n, • v + ) - v_(n 5 • v~) + ( u - “ pv + )S(<*>*)-^- J 
= ii 5 x (pXPr g e + • n 5 - Fr/e_ • n 5 ) 
and conservation of heat flux 

h s • (cpXV©| z= 0 * - V0| z=o -) = hs ■ [(cpt> + - v-)0| z =o + c(o u pv + - «-)] 

3<t> f 
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where c = c/( 1 - <r„), e is the rate-of-strain tensor (e± = e| z =o±), y is the surface-tension 
coefficient, tr u is the unbumt-to-bumt temperature ratio, and 5(<1> S ) and the unit norma n s are 
defined as 

5(<D,) = [l + O^j/S-c) 2 + (9 < l > i/ 3 )’) 2 ] 1/2 (13) 

n s = (-3<J>*/3jc, -d<t> s /3y, 1 )S(<P s ). 

Here, the factor multiplying y in equation (10) is the curvature —V • h., of the liquid/gas 
interface in the moving coordinate system, and the corresponding expressions for the gradient 
operator V and the Laplacian V 2 in this system are given by 
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However, the vector v still denotes the velocity with respect to the (x u x 2 ,x i ) coordinate 
system. 

We observe that the thermal and hydrodynamic fields are coupled through the temperature 
dependence of the mass burning rate A appearing in equation (9). When A is assumed to 
depend on pressure only, the strictly hydrodynamic problem for p, v and 4>< can be analyse 
separately [13, 14], In the present work, we wish to focus on this more extended coupling 
to determine how the hydrodynamic stability boundaries are modified when the local burning 
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rate depends on temperature as well. In connection with this we remark that although the mass 
burning rate A may be typically decomposed as 


A(p\ z =o+, 0|z=o) - Mp\z= o*. ®l 2 =o)exp 


~N( 1 

- a u 
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where A is a rate coefficient and N = EjR Q f a is the non-dimensional activation energy 
( E ) t this more explicit representation will not be needed in the stability analysis that follows. 
Instead, our (linear) stability results may be expressed in terms of the pressure and temperature 
sensitivities, defined as 
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respectively. However, we note that if the non-dimensional activation energy is large, then the 
last of equations (17) implies that Aq would likely be larger in magnitude than A p , a fact that 
will play a role in the relative scalings introduced for A p and Aq below. 

A non-trivial basic solution to the above problem, corresponding to the special case of a 
steady, planar deflagration, is given by 
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The remainder of the paper is devoted to a linear stability analysis of this solution. 


(18) 


3. The complete linear stability problem 


Prior to introducing any further approximations, we determine, in a standard fashion, the linear 
stability problem and the corresponding dispersion relation (in the form of a system of algebraic 
equations) that governs the behaviour of harmonic perturbations about the basic solution (18). 
In particular, in terms of the perturbation quantities (f> s (x y y,t ) = y> z, 0 — ^(Oi 
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0 _ ©o ( z )-<p s d 0 °/dz, the problem obtained when equations (l)-(3), ( 6 )-( 12 ) are linearized 
about the basic solution (18) is given by 
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where c = c( 1 — a u )~ 1 . 

Non-trivial harmonic solutions for <p s , u and C , proportional to e la " +lt| * +1 * 2V , that satisfy 
equations ( 1 9)— <2 1 ) and the boundary /boundedness conditions at z = ±00 are given by 
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where we have normalized the above solution by setting the coefficient of the harmonic 
dependence of <j> s to unity. Here, the signs of k\ and k 2 may be either positive or negative 

and we have employed the definition k — Jk\ + k\, and p, q , r and s are defined as 
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where we have used the fact, noted below equation (5), that pXPr g = pPri. 

Substituting this solution into the interface conditions (23)-(30) and using equation (19) 
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(dispersion relation) i co(k). In particular, these conditions are given by 
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While the above problem is linear in the coefficients b\-b%, the relationship for \a> is highly 
nonlinear. Accordingly, we seek asymptotic solutions for the neutral stability boundaries in 
the limit that the gas-to-liquid density and viscosity ratios p and /z are small, as is Fr 1 in 
the case of reduced gravity. We shall also introduce appropriate scalings for the sensitivity 
parameters A p and A©, where the limiting results for A© = 0 were derived previously [ 1 3, 14]. 
In particular, we introduce a bookkeeping parameter e <5C 1 and consider the realistic parameter 
regime 


p = p*e (i = Fr 1 


g ^ 0(1) 
g*€ - 0(0 


(52) 


where the upper (lower) choice for Fr~* corresponds to normal (reduced) gravity. We note 
that based on characteristic densities for liquid propellants and product gases at elevated 
temperatures, a realistic value for € is on the order of 10 -3 or 10 -4 . In this parameter regime, 
it was seen that the appropriate scaling for A p to describe the neutral stability region was 
A p = A* e, whereas the appropriate scale that describes the fully developed effects of thermal 
coupling turns out to be A© = Aj^e 1 ^ 4 , as shown in section 6 below. Thus, we observe that 
a 0 /a„ ~ 0(€-V*) > 175, which is at least several times as large as typical values of the non- 
dimensional activation energy but still consistent with the last of equations (17). However, 
as shown in section 7, the first significant effects of thermal coupling are in fact felt on an 
intermediate scale A© ~ 0(e 1/2 ), in which case A©/A p ~ 0(e -l/2 ) > 30, which is a quite 
typical value for N. Thus, the parameter regime of interest for describing modifications to the 
hydrodynamic stability boundaries arising from temperature-sensitivity effects coincides with 
realistic values of the sensitivity coefficients. 


4. The inviscid limit 


Although both liquid and gas-phase viscous effects were shown to be comparable in general 
to those due to surface tension when thermal coupling was neglected [13, 14], the qualitative 
nature of the cellular boundary remains preserved in the zero-viscosity limit, as does the 
existence of the pulsating boundary. Accordingly, we shall, for simplicity, first consider 
the effects of thermal sensitivity on hydrodynamic stability in the inviscid limit, reserving 
consideration of thermal coupling in the more algebraically involved viscous case for a 
future study (part II). In particular, before introducing the scalings indicated above, we first 
set Pri = Pr g = 0. Consequently, since q — ► oo according to equation (38), we set 
b 3 = b$ = b-i = 0 in equations (33M35), and observing that equations (48) and (49) then 
reduce to the no-slip conditions (43) and (44), we thus eliminate equations (41), (48) and (49) 
from the above set. Using the fact that r — > — i cop, the solution for the remaining coefficients 
in terms of i co is given by 
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Substitution of these results into equation (51) then gives a single implicit equation for ioj 
given by 
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where b%(\co) and b\oM are given explicitly by equations (53) and (57), provided A© ^ 0. 
Equation (59) is thus the dispersion relation that determines the neutral stability boundaries 
in that case. In the limit A© 0, the dispersion relation is given by setting the right-hand 
side of equation (57) to zero, in which case equations (58) and (59), which are then decoupled 
from equations (53)— (57), determine bg and b\o* Alternatively, we may solve equation (59) 
for b\o(ico) and substitute this result into equation (57), which then becomes the dispersion 
relation for arbitrary A©. We remark that the special case A© = 0 yields an analytically 
tractable solution for the neutral stability boundaries [12], but since all other results to be 
presented are obtained for the parameter regime defined by equations (52), we shall restrict 
ourselves to this regime for this limiting case as well. 


5. Hydrodynamic stability boundaries in the limit A© = 0 


In the limit A© = 0, the results obtained from the dispersion relation (57) in the parameter 
regime (52) are as follows [13]. First (see figure 1), the pulsating hydrodynamic stability 
boundary is given, to leading order in €, by A* = — p* for all disturbance wavenumbers (the 
exact relation is given by A* = — p*/( 1 — £p*)), with instability occurring below this critical 
value. For the cellular stability boundary, the corresponding leading-order expression depends 
on the magnitude of the wavenumber k. In particular, there are three wavenumber scales to be 
considered; the 0(1), or outer, scale k y a far outer scale kf — k€, and an inner scale defined by 
either k{ = k/e if Fr ~ 1 ~ 0(1) or k t — k/€ 2 if Fr -1 ~ 0(e). In each of these regions, the 
cellular stability boundary, on which the complex frequency ia> is identically zero, is given by 

A* ^ A* io) \p* ( 60 ) 

a ; ~ a *;" ~ ip*(p>*/ - 1) < 61 ) 

and 


A 


* 

p 


i * (0 


P*(P*S - h)/2ki 
P*(P*g* ~ ki)/2ki 


Fr~ l ~ 0(1) 
Fr~' ~ 0(e) 


(62) 


respectively. Matching these solutions to one another, a uniformly valid composite expansion 
A* p ic) (k) in the region A* < 0 (the basic solution can be shown to be unstable for A* > 0) 
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Figure 1. Inviscid hydrodynamic stability boundaries in the limit of zero thermal coupling 

{A e = 0 ). 


may be constructed as 

4*(c) 

a p 


A* (i) + /T„ (0) + A* if) - lim A* {i) - lim A* if) 

p P P ir^n^ P y 


V 


+ \ep' 2 yk + 


ep^g/lk 

e 2 p* 2 g’/2k 


(63) 


ki~+ oo y kf—*Q 

-* 2 ‘-m Fr~ x ~ 0(1) 

Fr~ l - 0(6) 

where the definitions of and k f have been used to express the final result in terms of k. We 

note that this asymptotic result, exhibited in figure 1, is far simpler than the exact result [12], 
which, in our present notation, is given by 

_ , P*6( 1 ~ P*€)Fr~ ] +p*€yk 2 - (1 - p*e)k 

p P p* 2 6 2(3 _ p*e)Fr~ l + p* 2 e 2 yk 2 + (1 


^0 


= 0. (64) 


p*e)(2- p*e)k 

It is thus easily seen from figure 1 that there is a stable region between the pulsating and 
cellular stability boundaries for A@ = 0, and this result is preserved (in fact, enhanced) when 
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viscous effects, which play a stabilizing role with respect to both boundaries for sufficiently 
large wavenumbers, are included in the analysis [13, 14]. We observe that for the Landau 
(cellular) form of hydrodynamic instability, gravity and surface tension are stabilizing with 
respect to sufficiently small and sufficiently large wavenumber disturbances, respectively. The 
essential difference, as illustrated in figure 1, between the normal and reduced gravity limits 
in the realistic parameter regime considered here is that in the latter instance, gravity is only 
capable of stabilizing disturbances whose wavenumbers are very small, 0(e 2 ), whereas in 
the former case, gravity is sufficiently strong to stabilize disturbances whose wavenumbers 
are 0(6). As a consequence, hydrodynamic instability becomes a long-wave instability 
phenomenon in the reduced gravity regime considered here, since, from equation (63), the most 
unstable wavenumbers are 0(^/6), rather than 0(1), in that case. We remark that although 
positive values of the pressure sensitivity A p are hydrodynamically unstable, zero and negative 
values of A p over certain pressure ranges are, in fact, characteristic of the so-called ‘plateau’ 
and ‘mesa’ types of solid propellants (cf Timnat [17]), as well as for the HAN-based liquid 
propellants mentioned above [15]. 

6. Hydrodynamic stability boundaries for As > 0 

Retaining the parameter scalings introduced in section 3, we first consider the 0(1) 
wavenumber regime. For the case As = 0 considered previously, it turns out [14] that 
the appropriate expansion of the dispersion relation i co(k) is of the form 

\(x) ~ 6~ ly/2 (i&>o + 6 l ^ 4 ia;] + 6^ 2 io>2 + • • •) (65) 

where the leading-order term was suggested by the explicit results available for the inviscid 
case, and the expansion in quarter powers of e is suggested by the leading-order term in the 
resulting expansion for p given below. We now introduce the scaling A© = 6" A@, where the 
value of n ^ 0 is to be determined so as to give a non-trivial modification in the dispersion 
relation at the appropriate order, corresponding to the first effects of a non-zero temperature 
sensitivity on the neutral stability boundary for 0(1) wavenumbers. 

Substituting equation (65) and the various scalings into equations (37), (40) and (53), 
expanded forms for p, s and are obtained as 

P ~ P0<r 1/4 + P\ + p 2 € 1/4 + • ■ • 5 ~ 0(6 1/2 ) (66) 

b% ^ C °°2p*k‘ k € ' - ^ (iw o)(iwi)e _3/4 - 3-[(i^) + 2(iw 0 )(ift>2)]f -l/2 + --- (67) 

where 

/7 0 = (iaj 0 ) l/2 Pi = \[\CO\ / (\CD 0 ) V2 + l] 

P 2 = |(io>o) _l/2 [l +4 k 2 +4i&>2 - (ift>i) 2 /ift>o]- 
Consequently, from equations (57) and (67), b io has the expanded form 

‘“'•‘"'jHFi 

+ ^|-^[(i^ o) 2 P* - * 2 ] - m )* 1/4 

+ ^ — ^ 2 ) ~ (' w i) 2 — 2(i&>o)(ia>2)jf '^ 2 + • • • J 

~<-1C+C< 1/4 + Ce ,/2 + '-] 


(69) 
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which, when substituted into the dispersion relation (59) along with the other expansions 
deduced thus far, yields a sequence of equations for the determination of the coefficients in 
the expansion (65). Thus, at leading order, 0(e“ rt “ 1/4 ), equation (59) yields b\°Q = 0, which, 
from the implied definition of b\°Q in equation (69) above, gives 

a * p p*~ if (ky ° )2p * - * 2 ] = 0 < 70a > 

or 


(iwo) 2 = 


(£)Wo-)- 


(706) 


Consequently, ia>o is identically zero for A* = — p*/ 2, which corresponds to the leading- 
order hydrodynamic cellular boundary (60) in the 0(1) wavenumber regime when A@ — 0. 
Values of A* > —p*/2 imply cellular instability, since in that case there exists, according 
to equation (706), a root of the dispersion relation such that Re{ia>o} > 0, Im{icuo} = 0. 
On the other hand, the stability of the region A* p < —p*/2 is still indeterminate, since 
equation (706) indicates that Re{ia>o} = 0 there. Hence, the next non-trivial term in the 
expansion of ico is required to determine whether this region is stable or unstable, although the 
fact that Im{ia>o) ^ 0 implies that disturbances have a pulsating character for values of A* 
below the cellular stability boundary. 

At the next order, 0(e“' 1 ), in the analysis of the dispersion relation, the equation 6 |q = 0 
is obtained, whence the definition of b\ l Q in equation (69) implies that iaq = 0. Proceeding 
with the next higher-order equation obtained at 0(e _n+1 < /4 ) from the expanded version of 
equation (59), it can be seen that terms other than those proportional to A@ -1 will appear 
provided the choice n = | is made. In that case, collecting terms of 0(e°) and using the 
previous results b\ = 6 ( t o = 0 leads to the equation 


-Pobfo +2k - ~[(\w 0 ) 2 p* +k 2 ] = 0 (71a) 

which, upon using the definitions of po and 6^, may be solved for ia >2 as 
lc ° 2 = 2p*k ~~ + P* + ^@(ia>o)~ 3/2 ] 

= ~Ti ( A *p~ P*)[K + P' + P* V2 k- {/2 A%(2A* p + p*)“ 3/4 ] (716) 

where the expression (706) for ia>o has been used to obtain the final equality. It is readily seen 
that for Aq = 0, Refill = 0 for A* = — p* < —p*/ 2. That is, as described in the previous 
section, there exists, since Im{ia>} 0, a pulsating neutral stability boundary that lies below 
the cellular boundary such that the region between the two. namely — p*/2 < A* < — p*, is 
stable [Re{io> 2 } < 0], and the region below the pulsating boundary is unstable [Refio^J > 0]. 

For Aq > 0, the pulsating boundary is modified. In particular, in the region A* < — p*/2 
below the cellular boundary, the principal value of the complex factor in equation (716) may be 
written as (2 A* +p*) -3 ' /4 — [-(2 A* + p*)]" 3, 4 e -3m/4^ an( j t j lus ^ neu t ra j stability condition 
Refia^} = 0 is given by 


K + P* - K2p^ /2 r I/2 A^[~(2A; r + /)]“ 3/4 0. (72) 

Writing A* p = — p*+<3p*/2, where a represents the deviation, in units of p*/2, in the pulsating 
stability boundary from its value in the zero-thermal- 'Sensitivity limit A^ = 0, the condition 
(72) becomes a — y' 2 p* ' 4 h “ 1 ' A *,( I • - (!) = 1 1, ( m eq u i v alent 1 y, 

a ( 1 n ri — 4p * k * A ( ,j . 


(73a) 
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Equation (73a) is an implicit expression for the pulsating stability boundary a(k), where the 
value a = 1 corresponds to the cellular boundary and a = 0 is the pulsating boundary for 
Aq = 0. Alternatively, writing A* — (— p*/2)(\ + b), where b = 1 - a is the negative 
deviation, in units of p* / 2, from the cellular stability boundary A* = — equation (73a) 

can be rewritten in terms of b as 


P( 1 - b) A = 4 p— l k- 2 Ag 

where b = 1 (b = 0) corresponds to the pulsating (cellular) boundary in the limit A* & = 0. 
A plot of jt(fr) for several different values of Aq is shown in figure 2, which, when rotated by 
—90° so that the fc-axis is horizontal, is readily interpreted in the context of figure 1, the lines 
A* = -p*/2 and A* = — p* in that figure corresponding to b = 0 and 1, respectively, in 

figure 2. 

Referring to figure 2, it is seen that for Aq > 0, the pulsating boundary becomes C-shaped 
(in the rotated frame of reference), the upper branch approaching the cellular boundary b — 0 
as k oo, and the lower branch approaching the original ( A@ =0) pulsating boundary b — 1 . 
These same limits are approached for any fixed value of k as A @ — ► 0. The portion within the 
C-shaped curve is the stable region, and thus not only is steady, planar burning intrinsically 
unstable for sufficiently small wavenumbers, but, for finite values of k, any crossing of the C- 
shaped boundary from the stable to the unstable region corresponds to the onset of a pulsating 
instability. As Aq increases, the turning point of the C-shaped pulsating boundary, which 
occurs at b = b c = f (i.e. at A* = -5p*/l) corresponding to the critical wavenumber 
k = k c = j£(2\ /p*) 1/2 A@ 2 ^ 21.83 Ag/p"t\ shifts to larger values of k as A^ increases. 
On the other hand, as A@ becomes small, the turning point shifts to small values of k such that 
k c eventually leaves the 0(1) wavenumber region for which equations (73) are valid. Thus, 
as A@ becomes small, the original pulsating and cellular boundaries are recovered in the 0(1) 
wavenumber regime, but as Aq becomes large, the original cellular boundary lies within the 
unstable region for 0(1) wavenumbers, and the basic solution becomes intrinsically unstable 
to oscillatory disturbances t- 

A composite asymptotic solution for the neutral stability boundary in the regime A© ~ 
0(t^ 4 ) is thus obtained by matching the cellular and pulsating boundaries in the far outer 
wavenumber regime, where the former is given by equation (61) and the latter is given trivially 
by A* = — p*, with the appropriate solution branch of equation (72) in the 0(1) wavenumber 
region. In particular, we denote the two solution branches of equation (72), which correspond 
to the portions of figure 2 that lie to the left and to the right of the minimum at b = b c , 
by A* io u) {k) and A* (o /) (fc), where the superscript ‘o’ denotes, as before, the outer, or 0(1), 
wavenumber region and the superscripts ‘u’ and ‘1* denote the upper and lower (rotate figure 2 
by —90°) portions of the double-valued pulsating boundary A* (k). These branches are given 

implicitly by equation (13b), where, by definition, b = — (2A*/p* + 1). Along the upper 
branch. A* (o ’ w) — > — p * /2 (i.e. b 0) as k -+ oo, which can be matched with equation (61) 
since, from equation (61), A* p {f) — ► —p*l 2 as kf -> 0. Similarly, A* p (oJ) — ► —p (i.e. 
b i) as k -> oo, which clearly matches the pulsating boundary A* = - p * in the far outer 
wavenumber region. As a result, a leading-order composite stability boundary spanning both 

t These results correct an erroneous conclusion drawn in a previous study [12], where an algebraic error in a similar, 
but less formal, calculation incorrectly suggested that the pulsating boundary recedes (to more negative values of A*) 
as increases, leaving the cellular boundary unaffected. 
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Pulsating Stability Boundary 
A e ~ 0(e 1 ' ,4 ) 1 k ~ 0(1) 



the outer and far outer wavenumber regions is given by 


a;w 




A*/’ u \k)+ep* 2 yk/2 

A* p l °'‘\k) 


a ; 2 —5p*/i 
A ; < -5p*/7 


(74) 


where the second term in the top expression has been expressed in terms of the outer 
wavenumber variable k. This composite boundary is shown in figure 3. We note, based 
on the above construction, that the lower branch of equation (74) is a pulsating boundary for 
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a 'p 



Figure 3. Composite hydrodynamic stability boundary for Aq ~~ 0(e 1 ^ 4 ). — 0.5, p* — y — 

1,6“ 0.005. 


all wavenumbers, whereas the upper branch transitions form a pulsating boundary for 0(1) 
wavenumbers to a cellular boundary for 0(€~ ] ) wavenumbers. Indeed, from equation (10b), 
we observe that in the outer wavenumber regime, the size of the upper region of oscillatory 
instability, which is bounded below by the upper branch of the pulsating stability boundary 
and above by the region of non-oscillatory instability beyond the old cellular boundary 
A* ~~ —p*f 2 for A® = 0, shrinks to zero as k becomes large. 


7. Evolution of the pulsating boundary for A@ <C 1 

The nature of the evolution, as Aq decreases, of the pulsating stability boundary depicted in 
figure 2 to that shown in figure 1 for the case of A@ = 0 may be determined by analysing the 
dispersion relation for smaller order-of-magnitude wavenumbers and appropriately rescaled 
values of A 0 . In particular, the results (70 b) and (73b) suggest that if we consider k ^ 0(0, 
then A@ ^ 0(e a/2 ) and \co 0 ~~ 0(e a ), where the latter imply that A @ and \(o are 0(6 (2a+1)/4 ) 
and 0(e (2a-!)/2 ), respectively. A preliminary analysis then suggests that qualitative deviations 
from the C-shaped nature of the pulsating boundary exhibited in figure 2 occur for wavenumbers 
of order € 1/2 , which is another intermediate wavenumber scale that lies between the inner and 
outer wavenumber regions that were defined in the previous section. The above analysis is 
thus repeated for k = ke xn and A 0 = A©* 1/2 , the other parameters remaining unchanged 
from their previously assumed orders of magnitude. Consequently, in place of equation (65), 
a solution for the dispersion relation is now sought in the form 


\(D ~ \co{) + € 1//4 ia>i + € X ^ 2 \a)2 + * * ■ • 


(75) 
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Substituting equation (75) and the new scalings introduced above into equations (37), (40) 
and (53), the new expanded forms of p , 5 and b% are given by 

P ~ Po + PiC /4 + p 2 € l/2 + ■■■ s~0(e) (76) 

, (ift»o ) 2 p*+k 2 1 . „ . „ ... 1 r ~ „ , n 

h - — f ; - T(ituo)(ift>i)c 7 ^[(ia;^) + 2(i&io)(i«2)]e 0 + • • • (77) 

2 p*k k 2k 

where 


Po = + (iwo) 1/2 ] p\ ~ \\cb\{\ + 4id>o) 1/2 p 2 = |i<i>2 ( 1 + 4ia>o)" 

As a result, equations (57) and (67) give the corresponding expansion of b\o as 

*10 ~ A + %G«0>V - 3 1/2 + • • 

A© l P* 2& 2 J 


1/2 


1 


+ y { — -j[(i&>o) 2 p* - P] - ^r(icbo)(\cO])€ 


1/4 


A© l 2& 


-T- { (i&>o) 2 p* — 3p} — (io?i) 2 — 2(i^)(ia> 2 )j 


e J/2 + < 


l (0) . ^ (1) v 4 r ( 2) 1/2 

°10 + 6 + t 


y 10 


10 ' 


+ * • 


(78) 


(79) 

Substituting these expansions into the dispersion relation (59), one obtains from the leading 
0(^°) terms in that equation the result bf^ = 0, or, according to equation (79), 


(i*o) 2 = ^ j (2 A* + p*). (80) 

Comparing equation (80) with equation (706), it is thus seen that the same leading-order result 
is obtained for icbo(k) for A© ^ 0(6 1/2 ) as was obtained for icoo(k) for A© ~ 0(6 1/4 ). In 
particular, \cbo is real and positive for A* > — p*/ 2, and is pure imaginary for A* ^ — p*/2. 
Thus, the prediction of cellular instability is recovered in the former case, but stability in the 
latter region can only be ascertained by calculating the next non-trivial term in the expansion 
(75). 

Collecting terms of 0(6 1/4 ) in the expanded version of equation (59), it is concluded that 
b% ~ 0, which, from equation (79), implies that \d>\ = 0. Proceeding to the next order, we 
collect terms of 0(6 1/2 ) in equation (59), which gives rise to the non-trivial equation 


p^ + ak + p'b™) 


V KUO / 


(81a) 


where, from equation (77), 6g 0) = — (i<^)(\co 2 )/k. Upon substituting the expressions for b ^ 
and (io>o) 2 from equations (79) and (80), we obtain the expression for \d >2 as 


id > 2 ~ k(A 




+p 


A<$p* 2 / 1 + 2 io>o — y 1 + 4i&>o \ "1 

k(2A* p + p*) V 1-2 c - y/\ + 4 ia>o / J 


( 81 *) 


where i<*o = (k/p*)(2A* p + p*) x > 2 = {k / p*)[-(2A* p + p*)] l/2 e i,r/2 , the second equality 
denoting the principal root when A* p < -p*/ 2. The expressions (81*) and (71*) for the 
two growth-rate corrections ia>2 and i<x>2 in their respective wavenumber regimes collapse to 
the same result in the limit of zero thermal sensitivity, yielding the same pulsating stability 
boundary A* = -p* as before. However, for non-zero A©, the forms of these two results 
differ, leading to qualitative differences in the corresponding neutral stability boundaries. In 
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particular, it will be shown below that the pulsating neutral stability boundary obtained from 
equation (81 b) approaches that shown in figure 1 (i.e. A* = — p*) in the limit of zero thermal 

sensitivity, whereas it approaches the form shown in figure 2 as A© becomes large. That is, 
the evolution between these two limiting forms occurs on the 0(6 l/2 ) wavenumber scale for 
0(e 1/2 ) values of the thermal sensitivity parameter A©. 

Setting the real part of i coi to zero, the equation governing the (pulsating) neutral stability 
boundary A p (k) in the present wavenumber regime is determined from equation (81b) as 




1 + 2i<£>o “ >/ 1 + 4ii>o 

k(2A* p + p*) ™ l 1 -2 c-y/\ + 4ia)o 


A&p 


Re 


(82) 


which, since ia>o depends on A p and k, is a highly nonlinear equation for A* (k). Thus, in 

order to determine the nature of the stability boundary as a function of A©, it is useful to 
consider certain limiting cases. We first consider possible solutions of equation (82) in the 
large (scaled) wavenumber regime k 1. In that limit, we observe from equation (70b) that 
ia>o ~ 0 (k) » 1, and thus a first approximation to equation (81) in this regime is given by 


A*+p*~ ^A®^l_Re!-v& ~0 (83a) 

p k(2A* p + p*) 

or, using the expression for the principal root icbo given below equation (81b), 

A; + p* - iV2p*V 2 k- l/2 A e [-(2A' p + p*j\~ V * ~ 0 (83 b) 

which, since k~ [/ 2 A & = (£/€ 1/2 )~ ,/2 (A e /6 !/2 ) = * _1/2 (A©/e 1/4 ) = fc _I/2 A@, is identical 
to equation (72). According to the discussion of equation (72) above, the turning point in the 
corresponding neutral stability boundary exhibited in figure 2 occurs at the value k c oc Aq 2 . 
Thus, we conclude that equation (83b) is a valid approximation for equation (82) in the large 
(scaled) wavenumber regime k 1 provided A^ 1, where the latter is a necessary condition 

to ensure that k is large everywhere along the stability boundary described by equation (83b). 
For such values of A©, the pulsating stability boundary for large k approaches that illustrated 
in figure 2 (with k replaced by k and A* e replaced by A©). 

In addition to the boundary described by equation (83b) and figure 2 for large k and 
A©, there is another pulsating branch that is present for small wavenumbers, as suggested 
by the observation that there exists a solution of equation (82) such that A* -► — p* 
as k 0, irrespective of the magnitude of A©. To simplify the analysis and further 
discussion, we mainly consider a simplified form of equation (82) by restricting further 
consideration to the r easonable limit in which the heat-capacity ratio c is small. Since 
1 + 2ia>o — s/T+^iu>q = (1 — y/l + 4ii>o ) 2 /2, a first approximation to equation (82) in that 
limit is given by 


A* + p* + „ Re{i/l +4ia>o — 1 1 ~ 0 c->0 

p 2k(2A* p + p*) I J 


(84) 


which clearly gives the same limiting approximation (83a) for large k. To analyse equation (84) 
for arbitrary k , it is convenient to convert it to an alternative form as follows. Writing 1 + 4io>o 
in the polar form 1 + 4io>o = rcostf + isintf, where r 2 = 1 4 - 16o>q, cos# = \/r and 
sin # = 4o>o/r, equation (84) may be rewritten as 

2^ ( a; + p*) = Mefyitd + l^^ + i] _ i] 


(85) 
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{d) Pulsating Stability Boundary 

k A e = 1 



(b) Pulsating Stability Boundary 

k Ae = 2/3 



unstable l 

o ' 2 oa o!e o.0 1 ^ 


Figure 4. Pulsating hydrodynamic stability boundaries for k — k€ 1/2 ^ 0(* 1/2 ) and decreasing 
values of A© = A©* 1/2 ~ 0(e^ 2 ). Jhe broken curves denote one- and two-term asymptotic 
expansions of the upper boundary for A© » 1, k » 1. (a) A© = 1; (b) A© = \ \ (, c ) A© = 
(d)A e =^(e)A e = \\{f)A e =^ 


where, from equation (80), Wq = — (k/p*) 2 (2A* p +p*). Substituting this expression for co^ into 
equation (85) and rearranging/squaring the latter twice so as to remove the fractional powers, 
a polynomial equation for the inverse relation Jc(A * p ; A 0 ) is obtained as 

4p(2 a; + P*) 3 (a; + p*) 4 - $>k 2 p* 2 A@(2A* p + p') 2 (A’ p + p*) 3 

+*[V 4 a|(2a; + p*)(a; + p*) 2 + p* 6 a£] - p* 6 a^(a; + p*) = o (86a) 

or, introducing b = -(2 A*/p* + 1) as before, 


Pb\l - b) 4 +4k 2 A @ b 2 (\ 


b ) 3 + *[5A|£(1 - b ) 2 - 4Aq/p*] + 2Aq(1 - b) = 0. 


( 86 *) 


An even more compact version is obtained by defining k = b( 1 — b)k/A e, giving the cubic 
equation 

P + 4* 2 + (5 - a)k + 2 = 0 (86c) 

where or(£) = 4A@/[p*£(l — £) 2 ]. However, since it is the relation A* (£; A©), or equivalently, 
^ (ft; A©), that is ultimately desired, we shall mainly consider the form given by equation (86/?). 
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In the regime A© » 1 there exists a solution k ~ 0(A|) as suggested above. More 
precisely, this solution of equation (86 b) may be sought in the form 

k ~ A 2 @ {k 0 + Ae*i + Aq 2 £ 2 + • ' • ) (87) 

whence, upon substitution of this expansion into equation (86 b), the coefficients k, are 
determined recursively as 

k 0 = 2b~ i/2 (\ -b)- 2 /p* in k = -2b-\\-bT x .... (88) 

As expected, the leading-order approximation k ~ A%k 0 is identical in form to the result (73b) 
and the next term in the above expansion thus provides a negative correction to this result. On 
the other hand, for large A© there also exists a solution k « 1 of the form 

k ~ Aq 1 (ko + Ag 1 k\ + Aq fa + • • ■ )• 

Substituting the latter into equation (86b), the k, are also detemined recursively as 

k 0 = p*(l -b)/2 ^i=0 k 2 = 5p* 2 b(l -b) 3 /8 .... (90) 

where the last of these provides a positive correction to the leading-order approximation 
k ~ (2A©) -1 p*(l - b). Although the branches described by equations (87) and (89) were 
obtained for A© » 1, they remain reasonably accurate representations of the actual branches 
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(«) 


Pulsating Stability Boundary 
A* = 1/4 


if) 


Pulsating Stability Boundaiy 
= 1/16 



even for 0(1) values of this parameter, as shown in figures 4(a) and ( b ) (discussed immediately 
below). 

* decreases through 0( 1 ) values, the turning point in the stability boundary for large 

k occurs for smaller values of k , and eventually the expansion for large k ceases to be a valid 
approximation of equation (86/?). Simultaneously, the branch corresponding to small values of 
k reaches larger values of k such that the above expansion for small k ceases to be valid as well. 
These two branches eventually intersect, restoring a region of stability for a range of negative 
values of A* . This development is illustrated in figures 4 which were obtained directly 
from the exact roots of the cubic equation (86 b) with respect to k (for comparison, we also 
indicate the asymptotic representation for the upper branch corresponding to one- and two-term 
truncations of equation (87)). We observe that after the two pulsating branches intersect (as, 
for example, in figures 4 (c)-(e)), there exist two critical values of k corresponding to the onset 
of pulsating instability at corresponding critical values of b, or equivalently, A*. That is, in 
this range of temperature sensitivities, a pulsating hydrodynamic instability occurs for either 
sufficiently small or sufficiently large negative values of A* as either the left or right (upper or 

lower in the rotated frame of reference) stability boundary is crossed. As A© decreases further, 
the separated branches move further apart until ultimately, in the limit A© -► 0, we recover 
from equations (86) the two roots b = 0 and 1 corresponding to the cellular (A* = —p* /2) 
and pulsating (A* = -p*) stability boundaries in the present (k ~ 0(e 1/2 )) wavenumber 
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regime. Thus, as the temperature-sensitivity parameter A e decreases to small values on the 
scale A&/A p — 0(e~ 1/2 ) — 50, the region of hydrodynamic stability depicted in figure 1 is 
recovered. 


8. Summary 

The present work has continued the formal asymptotic treatment of hydrodynamic instability 
in liquid-propellant combustion based on a generalized Landau/Levich model that allows for 
a dynamic dependence of the burning rate on local pressure and temperature perturbations. 
The focus in the present study was on the effects of thermal coupling, represented by non- 
zero values of an appropriately defined temperature-sensitivity parameter A 0 , on the nature 
of a pulsating stability boundary that exists for negative values of a corresponding pressure- 
sensitivity parameter A p . Again exploiting the smallness of the gas-to-liquid density ratio 
p, which provides the underlying basis for the asymptotic treatment of the stability problem, 
it was shown that for sufficiently large values of A 0 , the pulsating boundary possesses a 
turning point that renders steady, planar burning intrinsically unstable for sufficiently small 
wavenumber perturbations. As this parameter decreases, the stable region re-emerges such that 
in the limit that A © vanishes, the original separated pulsating and cellular stability branches 
are recovered. The evolution of the pulsating boundary between these two structures was 
shown to occur on an intermediate wavenumber scale relative to the inner, outer and far outer 
wavenumber scales that emerge in the asymptotic limit of small density ratios in the absence 
of thermal coupling. Values of the ratio Aq/A p over which this evolution occurs were shown 
to be roughly of the order of the overall activation energy, suggesting that a pulsating, rather 
than the classical Landau (cellular), form of hydrodynamic instability may be the more likely 
manifestation of hydrodynamic instability in at least some types of liquid propellants. 
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